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ABSTRACT 


This  paper  Indicates  how  exponential  random  variables  can  be  efficiently 
used  In  a  variety  of  simulation  problems.  One  of  the  problems  Is  the  simu¬ 
lation  of  order  statistics  from  a  normal  population.  Me  discuss  the  general 
problem  of  simulating  order  statistics  and  then  consider  the  normal  case. 

start  by  showing  how  the  Von-Neumann  rejection  approach  via  the  exponen¬ 
tial  can  be  modified  to  become  an  efficient  algorithm  for  generating  a  nor-  _ 
mal  and  then  present  a  method  for  generating  normal  order  statistics.  -We- 
show  how  to  use  the  exponential  to  efficiently  simulate  random  permutations 
with  weights.  We- then  consider  the  problem  of  simulating  a  2-dlmenslonal 
Poisson  process  both  for  a  homogeneous  and  nonhomo geneous  Poisson  process. 


SIMULATION  USES  OF  THE  EXPONENTIAL  DISTRIBUTION 


Sheldon  M.  Ross  and  Zvl  Schechner 


1.  SIMULATING  ORDER  STATISTICS 

Consider  the  problem  of  simulating  the  order  statistics  ^(1)  ^  *(2)  ^ 

...  <  from  a  sample  of  size  n  from  the  continuous  distribution  F  . 

If  F  Is  Invertible,  this  can  be  accomplished  by  simulating  ^  ...  < 

U,  ,  the  order  statistics  from  a  sample  of  n  uniform  (0,1)  random  varl- 
(n) 

ables  and  then  setting  ~  ^  generate  the 

simulate  n  uniform  (0,1)  random  variables  U^^ . U^  and  then  order  them. 

(It  should  be  mentioned  at  this  point  that  the  time  necessary  for  sorting  the 
U^'s  Is  not  nearly  as  large  as  has  been  Indicated  In  the  simulation  liter¬ 
ature  (see,  for  Instance,  [5]).  The  reason  being  that  the  sorting  need  not  be 
done  via  a  general  purpose  procedure  (such  as  quicksort  which  requires  on 
the  order  of  n  log  n  comparisons)  but  rather  can  be  done  with  a  procedure 
that  utilizes  the  fact  that  the  values  to  be  sorted  are  generated  from  a 
uniform  distribution  (see  5.2.1  of  [2]).  This  leads  to  an  algorithm  requir¬ 
ing  on  the  order  of  n  comparisons.)  We  now  present  another  approach  for 
generating  f  •••.  . 

Let  Xj^,  ...,  be  Independent  exponential  random  variables  with 

rate  1  and  Interpret  X^  as  the  1th  Interarrlval  time  of  a  Poisson  process. 


and  so  Is  the  time  of  the  1th  event.  Now  apply 


Now  set  S.  =  'l  X  and  so  S  Is  the  time  of  the  1th  event.  Now  apply 
j=l  ^ 

the  well-known  result  that  given  S  . ,  ,  S  ,  ...,  S  are  distributed  as  the 

n+1  1  n 

ordered  values  of  a  set  of  n  uniform  random  variables.  Hence, 


^ - ,  ,..,g -  have  the  joint  distribution  of  ,  ...»  •  That  Is, 

^n+1  n+1 

we  have  the  following  algorithm: 
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2 


Step.  1:  Generate  n  +  1  independent  exponential  random  variables 


X  . .  . 
n+1 


Step  2:  Set  U 


Remarks : 


i  /n+1 

I  x/  Y  X,i=l,  ...,n. 


The  above  procedure  was  discussed  In  [4].  The  exponential  random  vari¬ 
ables  can  be  generated  either  by  using  being  uniform  on 
(0,1),  or  by  using  any  of  the  other  known  algorithms. 

Suppose  now,  as  in  [6],  that  of  the  n  order  statistics,  »  •••» 


V)  ““xly  ••••  “(l+k) 


To  simulate  this  set  of 


k  +  1  of  the  n  order  statistics,  start  by  simulating  X  ,  a  gamma  random 
variable  with  parameters  (i,l)  (either  by  summing  exponentials  when  i  is 
small  or  by  using  one  of  the  algorithms  that  is  more  efficient  than  this  when 
i  is  large) .  Interpret  X  as  the  time  of  the  ith  event  of  a  Poisson  process 
with  rate  1.  Now,  simulate  k  exponential  random  variables  with  rate  1 — call 


them  Xj^,  . . . ,  — and  set 


X  ,  j  =  1 . k  .  Now,  simulate 


Y  ,  a  gamma  random  variable  with  parameters  (n  +  1  -  i  -  k  ,  1)  ,  which  will 
represent  the  time  between  the  (i  +  k)th  and  the  (n  +  l)st  event.  Then 


^l+k  ^  ’  ^l+k  " . ^l+k  ^  ^ 


will  have  the  desired  distribution. 


That  is,  we  have  the  following  algorithm  for  simulating  . 

^(i+k)  ' 

Step  1;  Generate  k  +  2  independent  random  variables  X  ,  Y  ,  X^^,  ...,  Xj^ 
with  X  being  gamma  with  parameters  (i,l)  ,  Y  being  gamma  with 
parameters  (n+l-l-k,l)  ,  and  X^  being  exponential  with 


rate  1. 
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I 

I 

I 

I 


m 

I 

I 


I 


I 


Step  2;  Set  for  j  *  0,1,  ...,  k 


X+  1  X 

U  _ _ 

^1+j)  k 

X  +  I  X  +  Y 

Even  easier  than  simulating  uniform  order  statistics  is  simulating  ex¬ 
ponentially  distributed  ones.  Let  . . . ,  Y^  be  independent  exponential 

random  variables  with  rate  1  and  set 

Y 

(1)  n 

^(2)  °  ^(1) 

^(j)  “  ^(J-1)  n  -  J  +  1 
^(n)  “  ^(n-1)  \  • 

It  then  follows  from  the  lack  of  memory  property  of  exponentials  plus  the 
fact  that  the  minimum  of  independent  exponentials  is  exponentially  distrib¬ 
uted  with  rate  equal  to  the  sum  of  the  rates  of  the  exponentials  over  which 
we  are  minimizing  that  •••»  ^(n)  distributed  as  the  order  sta¬ 

tistics  of  a  set  of  n  Independent  exponentials  with  rate  1. 

-1  -Y  -1  -Y 

Since  F  (1  -  e  )  and  F  (e  )  both  have  distribution  F  when  Y 

is  exponential  with  rate  1,  it  follows  that  we  can  simulate  a  set  of  n 

order  statistics  from  F  by  first  simulating  the  exponential  order  statis- 

-1  -Y 

tics  Xqj,  ...,  and  then  either  setting  «  F  (1  -  e  (1))  , 

i  ”  1,  ...,  n  or  setting  X^^  ••  F  ^(e  ,  i  =  1,  ...,  n  . 
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9 


9 


P 


9 


4 


‘*X 

Hence,  if  F  (1  -  e  )  or  F  (e  )  Is  as  easy  or  almost  as  easy  to 
calculate  as  F  ^(x)  ,  it  is  more  efficient  to  first  simulate  exponential 
rather  than  uniform  order  statistics. 


I 
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2.  SIMULATING  NORMAL  RANDOM  VARIABLES 

Let  Z  denote  a  unit  normal  random  variable  and  set  X  =  jz]  .  The 
Von-Neumann  rejection  method  for  simulating  X  ,  via  the  exponential  dis¬ 
tribution,  leads  to  the  following  well-known  algorithm: 

(a)  Generate  independent  exponentials  with  rate  1,  and  Y2  . 

(b)  Set  X  =  Y^  if  '^2  —  ^^1  ~  1)^/2  •  Otherwise,  return  to  (a). 

To  improve  on  this  note  that,  by  the  lack  of  memory  property  of  the  exponen¬ 
tial,  it  follows  that  if  we  reject  in  Step  (b) — that  is,  if  Y^  <  (Yj^  -  1)^/2 
or,  equivalently,  if  Y^  >  1  +  — then  Y^^  -  1  -  /2Y2  is  exponential  with 

rate  1.  Hence,  we  can  use  this  as  one  of  the  exponential  values  in  the  next 

iteration.  In  fact,  even  if  we  accept  in  (b) ,  we  can  generate  an  exponential 

2 

(independent  of  X)  by  computing  Y^  -  (Y^^  -  1)  /2  which  will  be  exponential 
with  rate  1. 

Hence,  summing  up,  we  have  the  following  algorithm  which  generates  an 
exponential  with  rate  1  and  an  independent  unit  normal  random  variable: 

Step  1:  Generate  Y^^  ,  an  exponential  random  variable  with  rate  1. 

Step  2:  Generate  Y^  ,  an  exponential  with  rate  1. 

Step  3:  If  Y^  -  (Y^  -  1)^/2  >  0  ,  set  Y  =  Y^  -  (Y^^  -  1)^/2  and  go  to 

Step  4.  Otherwise,  reset  Y^^  to  equal  Y^^  -  1  -  /2^  and  return 
to  Step  2. 

Step  4:  Generate  a  random  number  U  and  set 


if  U  1/2 


-Y^  if  U  >  1/2 


The  random  variables  Z  and  Y  generated  by  the  above  are  Independent  with 

Z  being  normal  with  mean  0  and  variance  1  and  Y  being  exponential  with 

rate  1.  (If  we  want  the  normal  random  variable  to  have  mean  p  and  variance 
2 

a  ,  just  take  p  +  oz.) 

Remarks : 

(1)  As  Is  well  known,  the  above  requires  a  geometric  distributed  number  of 
Iterations  of  Step  2  with  mean  /2e/ir  S;i.32  . 

(2)  The  final  random  number  of  Step  4  need  not  be  separately  simulated  but 
rather  can  be  obtained  from  the  first  digit  of  any  random  number  used 
earlier.  That  Is,  suppose  we  generate  a  random  nunber  to  simulate  an 
exponential,  then  we  can  strip  off  the  Initial  digit  of  this  random 
number  and  just  use  the  remaining  digits  (with  the  decimal  point  moved 
one  step  to  the  right)  as  the  random  number.  If  this  Initial  digit  Is 
0,1, 2, 3  or  4  (or  0  if  the  computer  is  generating  binary  digits),  then 
we  take  the  sign  of  Z  to  be  positive  and  take  it  negative  otherwise. 

(3)  If  we  are  generating  a  sequence  of  unit  normal  random  variables,  then 
we  can  use  the  exponential  obtained  In  Step  4  as  the  Initial  exponential 
needed  In  Step  1  for  the  next  normal  to  be  generated.  Hence,  on  the 
average,  we  can  simulate  a  unit  normal  by  generating  1.32  exponentials 
and  computing  1.32  squares  and  .32  square  roots. 

2.1  Simulating  Normal  Order  Statistics 

Suppose  we  want  to  simulate  the  order  statistics  of  n  Independent 
unit  normal  random  variables.  To  do  so, we  will  generate  m  Independent  ex¬ 
ponentials  with  rate  1 — Wj^,W2,  ...,  — and  use  the  above  algorithm  to  gen¬ 

erate  a  binomial  distributed  number  of  normals.  However,  so  as  to  eliminate 


•9 

1 

7 

the  need 

to  order  these  resultant  normals,  we  shall  first  order  the  exponen- 

i 

tlals.  We  thus  have  the  following  approach: 

• 

Step  1: 

Generate  Independent  exponential  random  variables  with  rate  1, 

1 

W, . W  and  set 

1’  ’  m 

• 

' 

“(1)  -  "i'" 

P 

*^(3)  ■  ^(j-1)  ■^m-j+1  ’  . “• 

• 

r 

C 

Step  2: 

Set  k  -  0  ,  H  *  1  . 

Step  3: 

Generate  Y  ,  an  exponential  random  variable  with  rate  1. 

R 

Step  A: 

If  k  ^  m  stop;  otherwise,  reset  k  to  equal  k  +  1  . 

it 

Step  5: 

If  Y - ®  ,  reset  Y  to  equal 

i 

Y  2 

P 

i 

reset  1 

to  equal  f  +  1  and  go  to  Step  4.  Otherwise,  if  Y - ^ -  <  0  , 

P 

go  to  Step  3. 

The 

result  of  the  above  will  be  the  random  variables  X^j  ,  ...» 

R 

where  R 

*  1  -  1  Is  blnomlally  distributed  with  parameters  m  and  p  == 

• 

Since  the  above  Is  equivalent  to  using  the  rejection  method  on  the 

m  Independent  and  identically  distributed  (lid)  random  variables  W^, 

1 

,  It  follows  that  the  accepted  variables  are  also  lid.  Hence,  given  R  , 

X(i),  .. 

,  X^^^  are  distributed  as  the  order  statistics  of  a  set  of  R 

random  variables  having  the  distribution  of  |z|  where  Z  Is  a  unit  normal. 

1 

P 

■ 

» 


However,  as  R  need  not  equal  n  ,  we  are  not  finished.  If  R  >  n  ,  then 
we  need  randomly  choose  R  -  n  of  the  indices  1,  R  and  then  eliminate 

the  variables  having  these  indices.  For  instance,  if  n  =  5  and  R  =  7  and 

the  random  choice  of  2  of  the  first  7  Integers  is  1  and  4,  then  the  resultant 

5  order  statistics  are  ,  i  =  1,  ...,  5  where  =  X^2)  »  ^(2)  ° 

X(3)  ,  “  ^(5)  ’  ^(4)  *  ^(6)  ’  ^(5)  ~  ^(7)  ‘  best  way  to  randomly 

choose  R  -  n  of  the  indices  1,  ...,  R  depends  on  the  relative  size  of 
R  -  n  in  comparison  to  R  .  If  R  -  n  is  small  compared  to  R  ,  then  we 

can  generate  random  numbers  U  and  eliminate  the  index  [RU]  +  1  continuing 

until  R  -  n  distinct  integers  have  been  removed.  If  R  -  n  is  moderate  in 
comparison  to  R  ,  then  there  is  also  a  very  simple  algorithm  for  the  selec¬ 
tion  (see  Section  3.4.2  of  [1]). 

If  R  <  n  ,  then  we  can  reapply  the  above  steps  (with  a  different  value 
for  m)  and  then  merge  the  two  sets  of  order  statistics.  If  this  results  in 
greater  than  n  order  statistics,  then  we  use  the  above  technique  to  randomly 
eliminate  certain  of  them.  If  the  combined  set  is  still  of  size  less  than  n  , 
we  again  reapply  the  above  steps  (with  again  a  different  value  for  m)  until 
finally  we  have  at  least  n  order  statistics. 

The  problem  of  choosing  the  appropriate  value  of  m  when  k  order  sta¬ 
tistics  need  be  generated  remains.  As  E(R]  =  .7576  m  and  /Var  R  =  .4285  /m  , 
we  can  be  almost  certain  of  having  R  ^  k  if  we  take  m  to  be  approximately 
equal  to  ^  ^  +  Ak/ 3  .  However,  perhaps  a  value  around  k  (which  would 
give  one  roughly  a  fifty-fifty  chance  of  having  R  ^  n)  would  work  out  better. 
Numerical  experience  will  be  necessary  to  determine  the  best  m  when  the  order 
statistics  of  k  of  the  variables  ire  needed. 
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Once  the  above  is  accomplished,  we  will  have  the  order  statistics  of 
the  absolute  values  of  n  unit  normals.  To  obtain  the  order  statistics  ^ 

of  n  unit  normals,  generate  n  uniforms — one  for  each  order  statistics — 
and  break  the  list  of  order  statistics  of  the  absolute  normals  into  2 

parts — those  having  a  U-value  less  than  1/2  and  those  having  it  greater  ^ 

than  1/2 — keeping  the  same  relative  ordering  within  the  2  lists.  The  re¬ 
versed  order  of  the  second  list,  with  each  element  In  this  list  given  a 
negative  sign,  followed  by  the  first  list  now  yields  the  order  statistics  ^ 

of  n  unit  normals. 

Remark : 

As  In  Section  2,  the  final  n  uniforms  needed  to  determine  the  signs 
of  the  unit  normals  need  not  actually  be  generated  but  can  be  "stripped  off" 
from  the  random  numbers  used  earlier. 


M 
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SIMULATING  RANDOM  PERMUTATIONS  WITH  WEIGHTS 


Consider  an  urn  containing  n  balls — the  1th  of  which  weighs  w^  , 
i  =  1,  n  .  The  balls  are  sequentially  removed,  without  replacement, 

from  the  um  according  to  the  following  scheme.  At  each  withdrawal,  a  given 
remaining  ball  will  be  removed  with  probability  equal  to  its  weight  divided 
by  the  sum  of  the  weights  of  the  remaining  balls.  We  are  Interested  in 
simulating  the  random  permutation  I^  =  ^^1’  •••»  where  I^  is  the  index 

of  the  ith  ball  to  be  removed. 

The  following  algorithm  for  simulating  I  is  suggested: 

Step  1:  Simulate  independent  exponential  random  variables  ...»  each 

having  rate  1. 


Step 

_2: 

Set 

^i  = 

Ef/Wi  ,  i  =  1,  . . . ,  n  . 

Step 

_2; 

Let 

I  = 

(I^,  ....  I^)  where  1^ 

of 

...  x^  . 

That  the  above  algorithm  indeed  works  follows  from  the  lack  of  memory 
property  of  the  exponential  distribution  in  conjunction  with  the  fact  that 
if  ,  i  ^  1  ,  are  independent  exponentials  with  respective  rates  , 


.  "l 

i\ -""(vv  ••••%)!•  1^- 


Remark; 

The  above  can  also  be  used  in  the  problem  of  [5]  which  is  concerned  with 
generating  a  weighted  random  sample  of  size  k  of  the  n  balls. 
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4.  SIMULATING  A  CIRCULAR  REGION  OF  A  TWO-DIMENSIONAL  POISSON  PROCESS 

A  point  process  consisting  of  randomly  occurring  points  in  the  plane  is 
said  to  be  a  two-dimensional  Poisson  process  having  rate  A  if 

(a)  the  number  of  points  in  any  given  region  of  area  A  is  Poisson  dis¬ 
tributed  with  mean  AA  ;  and 

(b)  the  nuid>er  of  points  in  disjoint  regions  are  Independent. 


For  a  given  fixed  point  0  in  the  plane,  we  now  show  how  to  simulate 
events  occurring  according  to  a  two-dimensional  Poisson  process  with  rate 
A  in  a  circular  region  of  radius  r  centered  about  0  .  Let  ,  i  ^  1  , 
denote  the  distance  between  0  and  its  ith  nearest  Poisson  point,  and  let 
C(a)  denote  the  circle  of  radius  a  centered  at  0  .  Then 


p|7rR^  >  b|  =  P{no  points  in  C(/b/Ti) }  =  e  . 


Also, 


=  p|r2  >  V(b  +  Trr^)/TT  [  R^^  =  r| 

=  pjno  points  in  C^V(b  +  nr^j/irj  -  C(r)| 
In  fact ,  the  same  argument  can  be  repeated  to  obtain 


Proposition  1;  With  R^  =  0  , 


2  2 

hR^  -  ,  i  ^  1  ,  are  Independent  exponentials  with  rate  A 
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In  other  words,  the  amount  of  area  that  need  be  traversed  to  encompass 
a  Poisson  point  is  exponential  with  rate  A  .  Since,  by  S3rmmetry,  the 
respective  angles  of  the  Poisson  points  are  independent  and  uniformly  dis¬ 
tributed  over  (0,2tt)  ,  we  thus  have  the  following  algorithm  for  simulating 
the  Poisson  process  over  a  circular  region  of  radius  r  about  0  : 

Step  1:  Generate  Independent  exponentials  with  rate  1,  .  stopping 


(  X  +  . . .  +  X  2 
N  »  min  jn  :  - - -  >  r 


Step  2;  If  N  =  1  ,  stop,  there  are  no  points  in  C(r)  .  Otherwise,  for 
i  =  1,  . . . ,  N  -  1  ,  set 


+  . . .  +  X^)  /ATI 


Step  3;  Generate  Independent  uniform  (0,1)  random  variables  • 

Step  A;  Return  the  N  -  1  Poisson  points  in  C(r)  whose  polar  coordinates 


(R^,2itU^)  ,  i  =  1,  ...,  N  -  1 


The  above  algorithm  which  requires  on  average  1  +  Airr  exponentials 
and  an  equal  number  of  uniform  random  numbers  can  be  compared  with  the 
procedure  suggested  in  [3],  This  latter  procedure  simulates  points  in 
C(r)  by  first  simulating  N  ,  the  nunber  of  such  points  and  then  uses  the 
fact  that,  given  N  ,  the  points  are  uniformly  distributed  in  C(r)  .  Hence, 
the  procedure  of  [3]  requires  the  simulation  of  N  ,  a  Poisson  random  varl- 
able  with  mean  Aur  and  must  then  simulate  N  uniform  points  on  C(r)  , 


- ^ 
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2  2 

by  simulating  R  from  the  distribution  ^^(a)  *  a  /r  and  0  from  uniform 
(0,2ti)  and  must  then  sort  these  N  uniform  values  in  Increasing  order  of  R  . 

The  main  advantage  of  our  procedure  is  that  it  eliminates  the  need  to  sort. 

The  above  algorithm  can  be  thought  of  as  the  fanning  out  of  a  circle 
centered  at  0^  with  a  radius  that  expands  continuously  from  0  to  r  .  The  - 

successive  radii  at  which  Poisson  points  are  encountered  is  simulated  by 
noting  that  the  additional  area  necessary  to  encompass  a  Poisson  point  is 
always.  Independent  of  the  past,  exponential  with  rate  X  .  This  technique  ^ 

can  be  used  to  simulate  the  process  over  noncircular  regions.  For  Instance, 
consider  a  nonnegatlve  function  g(x)  and  suppose  we  are  Interested  in 
simulating  the  Poisson  process  in  the  region  between  the  x-axis  and  g  with 
X  going  from  0  to  T  (see  Figure  1) . 


FIGURE  1 


To  do  so,  we  can  start  at  the  left-hand  end  and  fan  vertically  to  the  right 

a 

by  considering  the  successive  areas  J  g(x)dx  .  Now  if  X,  <  X  <  ... 

0  1  2 

denote  the  successive  projections  of  the  Poisson  points  on  the  x-axis,  then 

^i 

analogous  to  Proposition  1,  it  will  follow  that  (with  X_  ■  0)  X  /  g(x)dx  , 
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i  ^  1  ,  will  be  independent  exponentials  with  rate  1.  Hence,  if  we  simulate 
^l’^2’  •••»  independent  exponentials  with  rate  1,  stopping  at 


N  =  min 


n  :  +  •••  ^  ^  J*  g(x)dx| 


and  determine  by 


^1 

X  j  g(x)dx  «  Ej 
0 

^2 

^  j  g(x)dx  *  Ej 


Vl 

/ 

S-2 


g{x)dx  -  E 


N-l 


and  if  we  now  simulate  Uj^,  — independent  uniform  (0,1)  random 

numbers,  then  as  the  projection  on  the  y-axis  of  the  Poisson  point  whose  x- 
coordinate  is  is  uniform  on  (0,g(X^))  ,  it  follows  that  the  simulated 

Poisson  points  in  the  interval  are  (Xj^,U^g(X^) )  ,  i«  1,  ...,  N  -  1  . 

Of  course,  the  above  technique  is  most  useful  when  g  is  regular  enough 
so  that  the  above  equations  can  be  solved  for  the  X^  .  For  instance,  if 
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g(x)  *  y  (and  so  the  region  of  Interest  Is  a  rectangle) ,  then 


E  +  ...  +  Ej 

\  ^  .1=1.  ....  N  -  1 


and  the  Poisson  points  are 


(X^.yU^)  .1=1.  ....  N  -  1 


f 


5.  SIMULATING  TWO-DIMENSIONAL  NONHOMOGENEOUS  POISSON  PROCESSES 

Consider  now  a  nonhomogeneous  Poisson  process  with  Intensity  function 
X(x,y)  and  suppose  we  are  Interested  In  simulating  this  process  over  a 
region  R  .  Let  X  be  such  that  X(x,y)  £  X  for  all  (x,y)  e  R  .  A 
thinning  algorithm  which  first  simulates  a  Poisson  process  having  rate  X 
over  R  and  then  accepts  the  resulting  Poisson  point  (x,y)  with  prob¬ 
ability  — was  recommended  In  [3].  However,  we  recommend  a  conditional 
approach  that  first  simulates  N  ,  a  Poisson  random  variable  with  mean 

I'l  -  //  X(x,y)dxdy — which  we  assume  Is  calculable  (and  which  represents 

(x,y)eR 

the  nuiil>er  of  points  In  R) — and  then  chooses  N  points  In  the  region  R  by 
simulating  from  the  density  f(x,y)  =  ^ simulate  from  this  density, 
an  acceptance  rejection  procedure  that  simulates  a  point  (X,Y)  uniformly  In 
the  region  and  then  accepts  It  with  probability  X(X,Y)/X  can  be  used.  That 
Is,  we  have  the  following  algorithm: 

Step  1:  Simulate  N  ,  a  Poisson  random  variable  with  mean  |Rj  . 

Step  2 1  Simulate  a  point  (X,Y)  uniformly  distributed  In  R  and  a  uniform 
(0,1)  variable  U  and  accept  (X,Y)  If 

^  X(X,Y_). 

^  -  X 

If  there  have  been  a  total  of  N  points  that  have  been  accepted,  stop — 
otherwise,  return  to  Step  2. 

Remarks ; 

(1)  Each  random  point  (X,Y)  will  be  accepted  with  probability  |R|/XA(R) 

where  A(R)  ■  //  dydx  Is  the  area  of  R  . 

R 


Hence,  the  mean  number  of 
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iterations  of  Step  2  needed  to  generate  N  accepted  points  is 
NXA(R)/jRj  .  As  N  has  mean  jRj  ,  the  above  thus  needs  on  average 
XA(R)  iterations  of  Step  2. 

(ii)  The  method  for  simulating  (X,Y)  uniformly  in  R  depends  on  the 
geometry  of  R  .  One  possibility  is  to  enclose  R  in  a  rectangle 
and  thus  randomly  choose  a  point  (X,Y)  in  this  rectangle  and  then 
accept  it  if  (X,Y)  e  R  . 

(iii)  The  fanning  out  technique  of  Section  3  can  still  be  employed  when 
X(x,y)  is  easily  Integrated  and  R  is  "nice."  For  instance, 
suppose  R  is  the  region  of  Figure  1.  Then  if  X^  <  X^  <  . . . 
denote  the  successive  projections  of  the  Poisson  points  on  the 

i  g(x) 

x-axis,  then  with  X_  *  0  ,  s  s  X(x,y)dydx  are  independent 
exponentials  with  rate  1.  Also,  given  a  Poisson  point  on  the  line 


X  *  X^  ,  the  y-coordinate  is  distributed  according  to  density 


X(X^,y) 


g(X^) 

/ 

0 


X(X^,y)dy  . 


9 
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